/**
 * @file   Surface.h
 * @author HirasawaYui <yui@Ubuntu18-04>
 * @date   Sat Dec 12 10:59:23 2020
 * 
 * @brief  用网格和数值解表现曲面。
 * 
 * 
 */

#ifndef _SURFACE_H_
#define _SURFACE_H_

#include <fstream>
#include <iostream>
#include <cmath>
#include <sstream>

#include "assert.h"
#include <string>
#include "include/tinyxml2.h"
#include "FEMSpace.h"
#include "Dofs.h"
using namespace tinyxml2;
typedef Eigen::VectorXd VectorXd;

#define TEMPLATE template<Dimension DIM>

TEMPLATE
class Surface
{
public:
    VectorXd _u;
    Mesh<DIM>* _mesh;
    Element<DIM> *_element;
    Surface(Mesh<DIM>* mesh, Element<DIM>* element,VectorXd u);
    virtual void WriteVTKData(const std::string ofilename);
    virtual void WriteVTKData3D(const std::string ofilename);
    void WriteVTUData(const std::string ofilename);
    void WriteVTUData3D(const std::string ofilename);
};

TEMPLATE
Surface<DIM>::Surface(Mesh<DIM>* mesh, Element<DIM>* element,VectorXd u)
{
    _u = u;
    _mesh = mesh;
    _element = element;
}

TEMPLATE
void Surface<DIM>::WriteVTKData(const std::string ofilename)
{
    std::cout << "Create vtk file..." << std::endl;
    std::ofstream os;
    os.open(ofilename + ".vtk");
    os << "# vtk DataFile Version 2.0\n";
    os << "VTK from Cpp\n";
    os << "ASCII\n";
    os << "DATASET STRUCTURED_GRID\n";
    os.precision(12);
    os.setf(std::ios::fixed, std::ios::floatfield);
    std::cout << "        Writing node data ..." ;
    os << "DIMENSIONS " << _mesh->n_dofs_per_row() << " " << _mesh->n_dofs_per_col() << " " << 1;
    std::cout <<"OK!" <<std::endl;
    os << std::endl;
    os << "POINTS " << _mesh->n_dofs() << " float\n";
    std::cout << "        Writing element data ..." ;
    for (int i = 0; i < _mesh->n_dofs(); i++)
    {
	Dofs<DIM> d = _mesh->DofsofIndex(i);
	os << d[0] << " " << d[1] << " " << 0 << std::endl;
    }
	std::cout << "OK!" << std::endl;
    os << "POINT_DATA " << _mesh->n_dofs() << std::endl;
    os << "SCALARS Value float\n";
    os << "LOOKUP_TABLE default\n";
    for (int i = 0; i < _mesh->n_dofs(); i++)
    {
	os << _u(i) << std::endl;
    }
    os.close();
    std::cout << "        create vtk file:" << ofilename + ".vtk" << " OK!" << std::endl;
}

TEMPLATE
void Surface<DIM>::WriteVTKData3D(const std::string ofilename)
{
    std::cout << "Create vtk file..." << std::endl;
    std::ofstream os;
    os.open(ofilename + "_3D.vtk");
    os << "# vtk DataFile Version 2.0\n";
    os << "VTK from Cpp\n";
    os << "ASCII\n";
    os << "DATASET STRUCTURED_GRID\n";
    os.precision(12);
    os.setf(std::ios::fixed, std::ios::floatfield);
    std::cout << "        Writing node data ..." ;
    os << "DIMENSIONS " << _mesh->n_dofs_per_row() << " " << _mesh->n_dofs_per_col() << " " << 1;
    std::cout <<"OK!" <<std::endl;
    os << std::endl;
    os << "POINTS " << _mesh->n_dofs() << " float\n";
    std::cout << "        Writing element data ..." ;
    for (int i = 0; i < _mesh->n_dofs(); i++)
    {
	Dofs<DIM> d = _mesh->DofsofIndex(i);
	os << d[0] << " " << d[1] << " " << _u(i) << std::endl;
    }
	std::cout << "OK!" << std::endl;
    os << "POINT_DATA " << _mesh->n_dofs() << std::endl;
    os << "SCALARS Value float\n";
    os << "LOOKUP_TABLE default\n";
    for (int i = 0; i < _mesh->n_dofs(); i++)
    {
	os << _u(i) << std::endl;
    }
    os.close();
    std::cout << "        create vtk file:" << ofilename + "_3D.vtk" << " OK!" << std::endl;
}

template<>
void Surface<2>::WriteVTUData(const std::string ofilename)
{
   ///构造 vtu 文件的 xml 模板
		XMLDocument* doc1 = new XMLDocument();
				
		XMLDeclaration * declaration = doc1->NewDeclaration();
		doc1->InsertEndChild(declaration);

/**     添加xml文件头的第二种方法
		const char* str = "<?xml version=\"1.0\" ?>";
		doc1->Parse( str );
*/
	
		XMLElement * element_VTKFile,
		           * element_Piece,
		           * element_Points,
		           * element_point_DataArray,
		           * element_Cells,
		           * element_Cells_DataArray_1,
		           * element_Cells_DataArray_2,
		           * element_Cells_DataArray_3,
		           * element_PointData,
		           * element_PointData_DataArray;
		           
		XMLNode * node_VTKFile, 
		        * node_UnstructuredGrid,
		        * node_piece,
		        * node_point,
		        * node_point_DataArray,
		        * node_Cells,
		        * node_Cells_DataArray_1,
		        * node_Cells_DataArray_2,
		        * node_Cells_DataArray_3,
		        * node_PointData,
		        * node_PointData_DataArray;

		
		element_VTKFile = doc1->NewElement("VTKFile");
		element_VTKFile->SetAttribute( "type", "UnstructuredGrid" );
		node_VTKFile = doc1->InsertEndChild( element_VTKFile );
		
		node_UnstructuredGrid = node_VTKFile->InsertEndChild( doc1->NewElement( "UnstructuredGrid" ) );
		element_Piece = doc1->NewElement( "Piece" );
		node_piece = node_UnstructuredGrid->InsertEndChild( element_Piece );

		element_Points = doc1->NewElement( "Points" );
		node_point = node_piece->InsertEndChild( element_Points );

		element_point_DataArray = doc1->NewElement( "DataArray" );
		element_point_DataArray->SetAttribute( "type", "Float32" );	
		element_point_DataArray->SetAttribute( "NumberOfComponents", 3 );	
		element_point_DataArray->SetAttribute( "format", "ascii" );	
		node_point_DataArray = node_point->InsertFirstChild( element_point_DataArray );
		
		element_Cells = doc1->NewElement( "Cells" );
		node_Cells = node_piece->InsertEndChild( element_Cells );

		element_Cells_DataArray_1 = doc1->NewElement( "DataArray" );
		element_Cells_DataArray_1->SetAttribute( "type", "Int32" );	
		element_Cells_DataArray_1->SetAttribute( "Name", "connectivity" );	
		element_Cells_DataArray_1->SetAttribute( "format", "ascii" );	
		node_Cells_DataArray_1 = node_Cells->InsertEndChild( element_Cells_DataArray_1 );

		element_Cells_DataArray_2 = doc1->NewElement( "DataArray" );
		element_Cells_DataArray_2->SetAttribute( "type", "Int32" );	
		element_Cells_DataArray_2->SetAttribute( "Name", "offsets" );	
		element_Cells_DataArray_2->SetAttribute( "format", "ascii" );	
		node_Cells_DataArray_2 = node_Cells->InsertEndChild( element_Cells_DataArray_2 );

		element_Cells_DataArray_3 = doc1->NewElement( "DataArray" );
		element_Cells_DataArray_3->SetAttribute( "type", "UInt8" );	
		element_Cells_DataArray_3->SetAttribute( "Name", "types" );	
		element_Cells_DataArray_3->SetAttribute( "format", "ascii" );	
		node_Cells_DataArray_3 = node_Cells->InsertEndChild( element_Cells_DataArray_3 );
		element_PointData = doc1->NewElement( "PointData" );
		element_PointData->SetAttribute( "Scalars", "u_h" );
		node_PointData = node_piece->InsertEndChild( element_PointData );
		element_PointData_DataArray = doc1->NewElement( "DataArray" );
		element_PointData_DataArray->SetAttribute( "type", "Float32" );	
		element_PointData_DataArray->SetAttribute( "Name", "u_h" );	
		element_PointData_DataArray->SetAttribute( "format", "ascii" );
		node_PointData_DataArray = node_PointData->InsertEndChild( element_PointData_DataArray );
   /// vtu 文件模板构造完成，开始写入相关数据     

 		const std::string& vtufilename(ofilename);
   		
	    std::ostringstream s1;
		std::string str1,str2;
	    s1.precision(15);
	    s1.setf(std::ios::scientific, std::ios::floatfield);

   		std::cout << "Creat vtu file..." << std::endl;	
   		std::cout << "\tWriting node data ..." << std::flush;
   
        int n_node = _mesh ->n_dofs();

        int n_element = _mesh ->n_element();
		element_Piece->SetAttribute( "NumberOfPoints", n_node );	
		element_Piece->SetAttribute( "NumberOfCells", n_element );	

	    for (int i = 0; i < n_node; i++)
	    {   
            Dofs<2> d = _mesh->DofsofIndex(i);
		    s1.str("");
		    s1  << std::endl << d[0] << "  " << d[1] << "  " << "0.0";
		    str1 = s1.str();
 			node_point_DataArray->InsertEndChild( doc1->NewText( str1.c_str() ));
        }
        std::cout << " OK!" << std::endl;
   		std::cout << "\tWriting element data ..." << std::flush;

        int n_data = n_node;
        std::vector<int> count(n_data, 0);
        std::vector<double> val(n_data, 0);

		int skip = 0;
	    for (int i = 0 ; i < n_element; i++)
   	    {
        const int n_dofs = _element -> n_Dofs();
        std::vector<int> Idx = _mesh->NodeofEle(i);
        std::vector<Dofs<2> > temp(n_dofs);
        for(int i = 0;i < n_dofs; i++)
	        temp[i] = _mesh->DofsofIndex(Idx[i]);
        _element->SetDofsList(temp);
           switch (n_dofs)
           {
              case 3: 
 		         s1.str("");
		         s1  << std::endl
		             << _mesh -> NodeofEle(i)[0] << "  "
		             << _mesh -> NodeofEle(i)[1] << "  " 
		             << _mesh -> NodeofEle(i)[2] << std::endl;
		         str1 = s1.str();
			     node_Cells_DataArray_1->InsertEndChild( doc1->NewText( str1.c_str() ));

 		         skip = skip + 3;
 		         s1.str("");
		         s1  << std::endl << skip;
		         str1 = s1.str();
			     node_Cells_DataArray_2->InsertEndChild( doc1->NewText( str1.c_str() ));
                 node_Cells_DataArray_3->InsertEndChild( doc1->NewText( "\n5\n" ));
                 break;

              case 4:
 		         s1.str("");
		         s1  << std::endl
		             << _mesh -> NodeofEle(i)[0] << "  "
		             << _mesh -> NodeofEle(i)[1] << "  " 
		             << _mesh -> NodeofEle(i)[2] << "  " 
		             << _mesh -> NodeofEle(i)[3] << std::endl;
		         str1 = s1.str();
			     node_Cells_DataArray_1->InsertEndChild( doc1->NewText( str1.c_str() ));
 		         skip = skip + 4;
 		         s1.str("");
		         s1  << std::endl << skip;
		         str1 = s1.str();
			     node_Cells_DataArray_2->InsertEndChild( doc1->NewText( str1.c_str() ));
                 node_Cells_DataArray_3->InsertEndChild( doc1->NewText( "\n9\n" ));
                 break;

              default:
                 //Assert(false, ExcInternalError());
                 break;
           }
            for (int j = 0;j < _element ->n_Dofs();j ++)
    	    {
		        int index = _mesh -> NodeofEle(i)[j];
		        count[index] += 1;
		        //val[index] += u_h.value(mesh.point(k), *the_element);
                val[index] += _u[index];
    	    }
        }
        for (int i = 0;i < n_node;i ++)
           val[i] /= count[i];
   		for (int i = 0; i < n_node; i++)
   		{
 		    s1.str("");
		    s1  << std::endl << val[i];
		    str1 = s1.str();
		    node_PointData_DataArray->InsertEndChild( doc1->NewText( str1.c_str() ));
        }
		doc1->SaveFile( (vtufilename + ".vtu").c_str() );
		std::cout << " OK!" << std::endl;
		std::cout << "\tcreat vtu file: " << (vtufilename + ".vtu").c_str() << " OK!" << std::endl;
        
	    delete doc1;
}
template<>
void Surface<2>::WriteVTUData3D(const std::string ofilename)
{
		///构造 vtu 文件的 xml 模板
		XMLDocument* doc1 = new XMLDocument();
				
		XMLDeclaration * declaration = doc1->NewDeclaration();
		doc1->InsertEndChild(declaration);
		XMLElement * element_VTKFile,
		           * element_Piece,
		           * element_Points,
		           * element_point_DataArray,
		           * element_Cells,
		           * element_Cells_DataArray_1,
		           * element_Cells_DataArray_2,
		           * element_Cells_DataArray_3,
		           * element_PointData,
		           * element_PointData_DataArray;
		           
		XMLNode * node_VTKFile, 
		        * node_UnstructuredGrid,
		        * node_piece,
		        * node_point,
		        * node_point_DataArray,
		        * node_Cells,
		        * node_Cells_DataArray_1,
		        * node_Cells_DataArray_2,
		        * node_Cells_DataArray_3,
		        * node_PointData,
		        * node_PointData_DataArray;

		
		element_VTKFile = doc1->NewElement("VTKFile");
		element_VTKFile->SetAttribute( "type", "UnstructuredGrid" );
		node_VTKFile = doc1->InsertEndChild( element_VTKFile );
		
		node_UnstructuredGrid = node_VTKFile->InsertEndChild( doc1->NewElement( "UnstructuredGrid" ) );
		element_Piece = doc1->NewElement( "Piece" );
		node_piece = node_UnstructuredGrid->InsertEndChild( element_Piece );

		element_Points = doc1->NewElement( "Points" );
		node_point = node_piece->InsertEndChild( element_Points );

		element_point_DataArray = doc1->NewElement( "DataArray" );
		element_point_DataArray->SetAttribute( "type", "Float32" );	
		element_point_DataArray->SetAttribute( "NumberOfComponents", 3 );	
		element_point_DataArray->SetAttribute( "format", "ascii" );	
		node_point_DataArray = node_point->InsertFirstChild( element_point_DataArray );
		
		element_Cells = doc1->NewElement( "Cells" );
		node_Cells = node_piece->InsertEndChild( element_Cells );

		element_Cells_DataArray_1 = doc1->NewElement( "DataArray" );
		element_Cells_DataArray_1->SetAttribute( "type", "Int32" );	
		element_Cells_DataArray_1->SetAttribute( "Name", "connectivity" );	
		element_Cells_DataArray_1->SetAttribute( "format", "ascii" );	
		node_Cells_DataArray_1 = node_Cells->InsertEndChild( element_Cells_DataArray_1 );

		element_Cells_DataArray_2 = doc1->NewElement( "DataArray" );
		element_Cells_DataArray_2->SetAttribute( "type", "Int32" );	
		element_Cells_DataArray_2->SetAttribute( "Name", "offsets" );	
		element_Cells_DataArray_2->SetAttribute( "format", "ascii" );	
		node_Cells_DataArray_2 = node_Cells->InsertEndChild( element_Cells_DataArray_2 );
		element_Cells_DataArray_3 = doc1->NewElement( "DataArray" );
		element_Cells_DataArray_3->SetAttribute( "type", "UInt8" );	
		element_Cells_DataArray_3->SetAttribute( "Name", "types" );	
		element_Cells_DataArray_3->SetAttribute( "format", "ascii" );	
		node_Cells_DataArray_3 = node_Cells->InsertEndChild( element_Cells_DataArray_3 );
		element_PointData = doc1->NewElement( "PointData" );
		element_PointData->SetAttribute( "Scalars", "u_h" );
		node_PointData = node_piece->InsertEndChild( element_PointData );
		element_PointData_DataArray = doc1->NewElement( "DataArray" );
		element_PointData_DataArray->SetAttribute( "type", "Float32" );	
		element_PointData_DataArray->SetAttribute( "Name", "u_h" );	
		element_PointData_DataArray->SetAttribute( "format", "ascii" );
		node_PointData_DataArray = node_PointData->InsertEndChild( element_PointData_DataArray );
   /// vtu 文件模板构造完成，开始写入相关数据     
 		const std::string& vtufilename(ofilename);
   		int n_node = 0.0;
		int n_element = _mesh ->n_element();
		double val_1 = 0.0;
        double pointx,pointy;
        int kk = 0;
	    std::ostringstream s1;
	    std::string str1,str2;
	    s1.precision(15);
	    s1.setf(std::ios::scientific, std::ios::floatfield);
   		std::cout << "Creat vtu file..." << std::endl;	
		for (int i = 0 ; i < n_element; i++)
        {
          const int n_dofs = _element -> n_Dofs();
		  std::vector<int> Idx = _mesh->NodeofEle(i);
          std::vector<Dofs<2> > temp(n_dofs);
          for(int i = 0;i < n_dofs; i++)
	        	temp[i] = _mesh->DofsofIndex(Idx[i]);
          _element->SetDofsList(temp);
          switch (n_dofs) 
          {
            case 3: 
            	n_node += 3;
            	break;
            case 4: 
            	n_node += 4; 
            	break;
          }
        }    	
		element_Piece->SetAttribute( "NumberOfPoints", n_node );	
		element_Piece->SetAttribute( "NumberOfCells", n_element );	
   		std::cout << "\tWriting element data ..." << std::flush;
        int skip = 0;		
        
	    for (int i = 0 ; i < n_element; i++)
   	    {
           	const int n_dofs = _element -> n_Dofs();
        	std::vector<int> Idx = _mesh->NodeofEle(i);
        	std::vector<Dofs<2> > temp(n_dofs);
        	for(int i = 0;i < n_dofs; i++)
	        	temp[i] = _mesh->DofsofIndex(Idx[i]);
        	_element->SetDofsList(temp);
           switch (n_dofs)
           {
              case 3: 
                 for (int k = 0; k < 3; k++)
                 {
                    kk = _mesh -> NodeofEle(i)[k];
					Dofs<2> t = _mesh ->DofsofIndex(kk);
                    pointx = t[0];
                    pointy = t[1];
                    val_1 = _u[kk];                                          
 		            s1.str("");
		            s1  << std::endl << pointx << "  " << pointy << "  " << val_1;
		            str1 = s1.str();
 			        node_point_DataArray->InsertEndChild( doc1->NewText( str1.c_str() ));
 		            s1.str("");
		            s1  << std::endl << val_1;
		            str1 = s1.str();
		            node_PointData_DataArray->InsertEndChild( doc1->NewText( str1.c_str() ));
                 }
		         s1.str("");
		         s1  << std::endl
		             << skip + 0 << "  "
		             << skip + 1 << "  " 
		             << skip + 2 << std::endl;
		         str1 = s1.str();
			     node_Cells_DataArray_1->InsertEndChild( doc1->NewText( str1.c_str() ));
				 skip = skip + 3;
 		         s1.str("");
		         s1  << std::endl << skip;
		         str1 = s1.str();
			     node_Cells_DataArray_2->InsertEndChild( doc1->NewText( str1.c_str() ));
                 node_Cells_DataArray_3->InsertEndChild( doc1->NewText( "\n5\n" ));
                 break;
              case 4:
                 for (int k = 0; k < 4; k++)
                 {
                    kk = _mesh -> NodeofEle(i)[k];
					Dofs<2> t = _mesh ->DofsofIndex(kk);
                    pointx = t[0];
                    pointy = t[1];
                    val_1 = _u[kk];                                          
 		            s1.str("");
		            s1  << std::endl << pointx << "  " << pointy << "  " << val_1;
		            str1 = s1.str();
 			        node_point_DataArray->InsertEndChild( doc1->NewText( str1.c_str() ));
 		            s1.str("");
		            s1  << std::endl << val_1;
		            str1 = s1.str();
		            node_PointData_DataArray->InsertEndChild( doc1->NewText( str1.c_str() ));
                 }
		         s1.str("");
		        s1  << std::endl
		             << skip + 0 << "  "
		             << skip + 1 << "  " 
		             << skip + 2 << "  "
		             << skip + 3 << std::endl;
		         str1 = s1.str();
			     node_Cells_DataArray_1->InsertEndChild( doc1->NewText( str1.c_str() ));

				 skip = skip + 4;
 		         s1.str("");
		         s1  << std::endl << skip;
		         str1 = s1.str();
			     node_Cells_DataArray_2->InsertEndChild( doc1->NewText( str1.c_str() ));
                 node_Cells_DataArray_3->InsertEndChild( doc1->NewText( "\n9\n" ));
                 break;
 

              default:
                break;
           }
        }
        std::cout << " OK!" << std::endl;
	doc1->SaveFile( (vtufilename + "_3D.vtu").c_str() );
	std::cout << "\tcreat vtu file: " << (vtufilename + "_3D.vtu").c_str() << " OK!" << std::endl;
        
	delete doc1;
}
#undef TEMPLATE
int creatpvd( std::vector< double > time_t, 
			  std::vector< std::string > filename,
			  std::string filenamepvd)
{
	
	if( time_t.size() != filename.size() )
	{
		std::cout << "number of time step is not equal to number of data files"  << std::endl;
		return 1;
	}
	XMLDocument* doc = new XMLDocument();
			
	XMLDeclaration * declaration = doc->NewDeclaration();
	doc->InsertEndChild(declaration);

		XMLElement * element_VTKFile,
		           * element_Collection,
		           * element_DataSet;
		           
		XMLNode * node_VTKFile, 
		        * node_Collection;

		element_VTKFile = doc->NewElement("VTKFile");
		element_VTKFile->SetAttribute( "type", "Collection" );
		node_VTKFile = doc->InsertEndChild( element_VTKFile );

		element_Collection = doc->NewElement( "Collection" );
		node_Collection = node_VTKFile->InsertEndChild( element_Collection );

	for(int i = 0; i < time_t.capacity(); i++ )
	{
		element_DataSet = doc->NewElement( "DataSet" );
		element_DataSet->SetAttribute( "timestep", time_t.at(i) );	
		element_DataSet->SetAttribute( "group", "" );	
		element_DataSet->SetAttribute( "part", 0 );	
		element_DataSet->SetAttribute( "file", (filename.at(i)).c_str() );	
		node_Collection->InsertEndChild( element_DataSet );
	}
	doc->SaveFile( (filenamepvd + ".pvd").c_str() );
	std::cout << "\t creat pvd file: " << (filenamepvd  + ".pvd").c_str() << " OK!" << std::endl;
	delete doc;

	return 0;

}

/**********************************************************/
/******************** EazyMesh Surf ***********************/
/**********************************************************/
class Surface_EM:public Surface<2>
{
public:
    Surface_EM(Mesh<2>* mesh, Element<2>* element,VectorXd u);
    void WriteVTKData(const std::string ofilename);
    //void WriteVTKData3D(const std::string ofilename);
};

Surface_EM::Surface_EM(Mesh<2>* mesh, Element<2>* element,VectorXd u):Surface<2>(mesh, element, u)
{
}   

void Surface_EM::WriteVTKData(const std::string ofilename)
{
    std::cout << "Create vtk file..." << std::endl;
    std::ofstream os;
    os.open(ofilename + ".vtk");
    os << "# vtk DataFile Version 2.0\n";
    os << "VTK from Cpp\n";
    os << "ASCII\n";
    os << "DATASET UNSTRUCTURED_GRID\n";
    os.precision(12);
    os.setf(std::ios::fixed, std::ios::floatfield);
    std::cout << "        Writing node data ..." ;
    os << std::endl;
    os << "POINTS " << _mesh->n_dofs() << " double\n";
    for (int i = 0; i < _mesh->n_dofs(); i++)
    {
	Dofs<2> d = _mesh->DofsofIndex(i);
	os << d[0] << " " << d[1] << " " << 0 << std::endl;
    }
    os << std::endl;
    os << "POINT_DATA " << _mesh->n_dofs() << std::endl;
    os << "SCALARS Value double\n";
    os << "LOOKUP_TABLE default\n";
    for (int i = 0; i < _mesh->n_dofs(); i++)
    {
	os << _u(i) << std::endl;
    }
    os << std::endl;
    std::cout << "OK!" << std::endl;
    std::cout << "        Writing element data ..." ;
    os << "CELLS " << _mesh->n_element() << " " << 4*_mesh->n_element() << std::endl;
    std::vector<int> tmp_ele;
    for (int i = 0; i < _mesh->n_element(); i++)
    {
	tmp_ele = _mesh->NodeofEle(i);
	os << 3 << " " << tmp_ele[0] << " " << tmp_ele[1] << " " << tmp_ele[2] << std::endl;
    }
    os << std::endl;
    os << "CELLS_TYPES " << _mesh->n_element() << std::endl;
    for (int i = 0; i < _mesh->n_element(); i++)
	os << 5 << std::endl;
    os.close();
    std::cout << "        create vtk file:" << ofilename + ".vtk" << " OK!" << std::endl;
    }

#else
#endif
